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ABSTRACT 


This  thesis  investigates  the  effect  of  the  ocean  water  attenuation  on  a  laser  beam 
fired  upward  inside  the  ocean.  The  laser  beam  spreading  due  to  scattering  is 
approximated.  The  method  used  is  a  computer  Monte  Carlo  simulation.  Angular 
spreading  of  light  caused  by  refraction  at  the  sea  surface  is  also  studied  and  compared 
with  the  ocean  results.  The  method  is  to  simulate  geometrical  light  rays  passing  through  a 
randomly  realized  ocean  surface  wave  model  and  derive  statistics  of  the  angular 
refraction.  The  results  of  this  work  can  be  used  for  detection  of  objects  in  the  water  and 
laser  communication  to  submerged  objects  from  an  airborne  or  space  platform. 
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I.  INTRODUCTION 


A,  OVERVIEW 

Laser  detection  of  objects  in  the  water  and  laser  communication  to  submerged 
objects  from  an  airborne  or  space  platform  requires  that  the  laser  beam  traverse  the  air- 
ocean  surface  on  the  way  to  the  object.  In  addition,  the  scattered  energy  of  a  source  below 
the  surface  must  traverse  the  ocean  on  the  upward  return  path.  The  rough  ocean  surface 
spreads  the  laser  energy.  Once  the  laser  beam  has  entered  the  ocean,  the  turbid  water 
scatters  the  energy,  further  spreading  the  laser  beam.  The  relative  significance  of  these 
two  processes,  and  the  angular  spread  of  a  laser  beam,  depends  on  the  depth  of  the  object, 
the  ocean  turbidity  and  the  size  of  the  beam  on  the  ocean  surface. 

B,  THESIS  OBJECTIVE 

To  investigate  the  beam  spread  versus  angle  of  incidence,  this  thesis  proposes  to 
simulate  the  angular,  refractive  effects  of  laser  energy  passing  through  the  ocean  surface 
and  the  scattering  within  the  ocean.  The  dependence  of  the  beam  size  on  the  ocean 
surface  and  depth  of  the  object  below  the  surface  will  be  considered.  The  primary 
approach  will  be  to  simulate  the  ocean  surface  with  different  models  of  the  spectrum  of 
the  ocean  surface  waves.  Analogous  to  atmospheric  wave-optical  laser  propagation 
codes,  a  representation  of  the  ocean  surface  can  be  simulated  by  starting  with  Gaussian 
random  numbers  low  pass  filtered  by  realistic  ocean  surface  wave  spectra.  Fourier 
transforming  back  into  x,  y  coordinates  provides  a  simulated  ocean  surface.  Using  ray- 
optics,  the  trajectories  of  a  sequence  of  optical  rays  will  reveal  the  angular  dispersion  of 
the  ocean  surface  waves.  Both  1-D  and  2-D  spectral  representations  of  the  ocean  waves 
will  be  considered.  After  the  refractive  effects  of  the  ocean  surface  waves  are  considered, 
the  effects  of  the  volume  scattering  will  be  included  using  the  known  scattering  functions 
for  the  ocean. 
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The  final  goal  is  to  develop  a  realistie  probability  distribution  for  the  angular 
spreading  of  the  laser  beam  including  the  turbid  scattering  as  a  function  of  the  depth  of 
the  object.  This  is  coupled  with  the  ocean  surface  refractivity  and  the  size  of  the  laser 
beam  footprint  on  the  ocean  surface. 

C.  THESIS  OUTLINE 

The  structure  of  this  thesis  is  organized  into  the  introduction,  background 
(Chapter  II)  and  four  additional  chapters.  Chapter  III  examines  the  wind-generated  ocean 
surface.  Starting  from  an  already  known  spectrum,  a  surface  model  suitable  for  the  scope 
of  this  thesis  will  be  approximated.  In  Chapter  IV,  the  air  water  interaction  at  the  ocean 
surface  will  be  examined,  using  Snell’s  Law.  Chapter  V  contains  the  analysis  of 
scattering  and  absorption  phenomena  inside  the  water  body.  A  Monte  Carlo  simulation 
will  be  used.  Next,  refraction  at  the  surface  will  be  again  evaluated,  this  time  including 
scattering  and  absorption.  Finally,  in  Chapter  VII  the  conclusions,  based  on  the  results 
obtained  from  the  analysis  in  the  previous  chapters,  are  presented. 
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II.  BACKGROUND 


In  this  chapter,  some  of  the  background  knowledge,  terminology  and  basie 
eoneepts  of  optical  oceanography  are  introduced. 

As  light  reaehes  the  earth’s  surfaee,  there  is  a  70%  probability  of  striking  a  water 
surfaee  and  interaeting  with  water.  At  this  interfaee,  various  phenomena  may  oeeur:  A 
brief  deseription  of  them  will  be  presented  next. 

A.  INHERENT  AND  APPARENT  OPTICAL  PROPERTIES 

Inherent  optieal  properties  are  those  whose  values  are  not  affeeted  by  the 
distribution  of  radianee  at  any  given  point.  Among  them,  the  most  important  eoneerning 
our  work  are  the  absorption  eoeffieient  (a),  the  scattering  eoefficient  (b)  and  the 

normalized  volume  scattering  function  ^(0) .  The  nature  of  the  light  field  inside  a  water 
mass,  given  the  incident  light  radiation,  is  a  function  of  these  properties. 

The  apparent  optieal  properties,  in  eontrast  with  the  inherent  ones,  are  those 
whose  value  at  any  given  point  in  the  medium  depends  on  the  radianee  distribution  at  that 
point.  These  eannot  be  regarded  as  properties  of  the  water  itself  but  of  the  light 
distribution  inside  the  water.  Typieally,  the  vertieal  attenuation  eoeffieient  is  examined 
for  the  downward  irradiance  and  the  irradianee  refleetanee. 

It  is  possible  to  ealeulate  the  inherent  properties  from  readily  observed  and 
measured  apparent  properties.  The  way  in  whieh  the  observed  properties  depend  on 
fluctuations  of  the  inherent  properties  ean  also  be  estimated. 

A  detailed  deseription  and  mathematieal  expression  of  all  the  above  properties 
and  the  fundamental  quantities,  whose  knowledge  is  essential  for  this  work,  is  presented 
in  Appendix  A. 

B,  REFRACTION  AT  THE  SEA  SURFACE 

The  interfaee  between  the  sea  water  and  the  air  ineludes  two  media  with  different 
optieal  densities.  Due  to  the  different  veloeity  of  propagation  in  these  two  media,  light 
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that  penetrates  in  a  water  body  changes  its  direction.  This  is  according  to  Snell’s  Law: 
sin  (j)^  =  sin  (j)^  ,  as  shown  in  Figure  1 . 


Figure  1 .  Snell's  Law  and  the  reflection  law 

where  cp^  and  (p^  correspond  to  the  angle  of  incidence  and  refraction,  respectively. 

and  are  the  refractive  indices  of  air  and  water  respectively.  These  indices  vary 

with  temperature,  salinity  and  wavelength.  For  all  practical  purposes,  the  above  effects 
can  be  neglected  by  assuming  that  water  has  a  constant  refractive  index  of  1.33.  The 
above  law  holds  for  both  directions  (light  entering  or  exiting  the  water  surface). 

C.  REFLECTION  AT  THE  SEA  SURFACE 

According  to  the  reflection  law,  the  light  beam  is  reflected  at  the  same  angle  with 
the  angle  of  incidence,  as  shown  in  Figure  1.  The  reflected  electric  vector  generally  has 
one  component  parallel  and  one  perpendicular  to  the  plane  of  incidence,  given  by 
Fresnel’s  equations. 
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Reflection  at  the  sea  surface  depends  on  the  solar  elevation,  the  existence  of 
waves  on  the  sea  surface  and  on  the  wavelength  (dispersion).  Furthermore,  in  bad 
weather  conditions,  it  is  difficult  to  discriminate  the  true  reflectance  from  the  short  wave 
back-scattered  light  from  the  sea  surface.  Figure  2  exhibits  the  solar  elevation 
dependence  of  the  reflection.  The  parameter  n  is  the  significant  ratio  of  sky  radiation  to 
global  radiation. 


Figure  2.  Amount  of  reflected  Energy  as  a  function  of  the  solar  elevation  for  various 

values  of  the  ratio  n  of  sky  radiance  to  global  radiation.  (From  [1]) 

1.  Sun  Glitter 

Glitter  is  a  special  reflection  phenomenon.  It  takes  place  when  a  flat  water  surface 
is  roughened  by  the  wind.  The  sun  image,  formed  by  specular  reflection,  has  many 
distinct  glittering  points.  As  the  water  surface  becomes  rougher,  the  phenomenon  is  more 
obvious  and  the  glittering  band  increases.  The  phenomenon  was  used  by  Cox  and  Munk 
for  estimating  the  sea  surface  slope. 
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2,  Dispersion  of  Reflection 

For  angles  lower  than  30  degrees,  observations  show  that  irradianee  refleetion 
depends  on  the  wavelength.  This  dispersion  is  oommon  to  other  phenomenon  eoneeming 
the  light  and  water  interaetion,  sueh  as  seattering  and  attenuation. 

D,  SCATTERING 

Seattering  and  absorption  are  the  two  main  proeesses  that  determine  light 
propagation  inside  the  oeean.  Seattering  involves  the  abrupt  deviation  of  light  from 
reetilinear  propagation.  It  ehanges  the  direetion  of  photon  transport,  “dispersing”  them  as 
they  penetrate  a  sample,  without  ehanging  their  wavelength.  The  seattering  eoeffieient, 
b(>.),  is  equal  to  the  fraetion  of  energy  dispersed  from  a  light  beam  per  unit  of  distanee 
traveled  in  a  seattering  medium,  in  m  ' . 

The  most  important  parameter  eoneeming  seattering  is  the  volume  seattering 
funetion,  whieh  represents  seattering  as  a  funetion  of  the  seattering  angle.  Important 
quantities  eoneeming  seattering  are  presented  in  Appendix  A. 
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Figure  3.  A  beam  of  light  with  Intensity  O.  enters  a  medium  of  volume  AV.  Part 

of  the  light  is  absorbed  inside  the  medium  (Oa),  part  of  the  light  is  scattered  at 
various  angles  'P  and  a  part  is  transmitted  outside  the  medium  (Ot).The  exiting 
rays  create  an  elliptical  footprint.  (After  [2]) 


1,  Geometric  Optics’  Scattering 

The  scattering  process  is  the  result  of  the  following  procedures;  Diffraction, 
where  light  is  reradiated  from  rectilinear  propagation;  Refraction,  where  light  penetrates 
the  object  and  emerges  with  or  without  internal  reflections;  Finally,  reflection,  where 
light  does  not  penetrate. 

2,  Scattering  in  Sea  Water 

Scattering  in  sea  water  is  produced  by  the  water  itself  and  the  suspended  particles 
as  well.  Since  small  particles  produce  wide  angle  scattering,  the  scattering  by  particles  is 
far  more  important.  Particle  scattering  is  mainly  dependent  on  the  particle  concentration, 
particle  size  and  optical  radiation  wave  length.  The  major  effect  comes  from  particles 
with  size  greater  than  the  wavelength  of  light.  Absorption  should  be  taken  into 
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consideration  as  well.  Diffraetion  is  largely  independent  of  the  eomposition  of  the 
partiele,  where  refraetion  and  refleetion  are  determined  by  the  refraetive  index  of  the 
partiele. 


3.  Scattering  by  Water  Background 

a.  Rayleigh  Theory 

Rayleigh  seattering  ean  be  eonsidered  as  a  moleeular  problem.  Aeeording 
to  Rayleigh  theory,  a  homogenous  eleetrieal  field  E  induees  in  eaeh  partiele  an  oseillating 
dipole,  whieh  radiates  in  all  direetions.  In  the  ease  of  N  partieles  with  a  smaller  size 
eompared  to  the  wavelength  of  the  incident  radiant  energy,  the  radiant  intensity  in  the  0 
direetion  (angle  respeetively  to  the  horizontal)  is  given  by  [3] 


2^2 


/  = 


8;z-"N«"E 


-(l  +  eos"6')  W/St 


(3.1) 


where  a  is  the  polarizability  of  the  partiele. 


b.  Fluctuation  Theory 

This  theory  relates  seattering  with  fluetuations  in  density  of  eoneentration 
in  a  small  volume  element  of  the  fluid.  Regarding  intensity  distribution  and  polarization, 
seattering  has  the  same  properties  with  Rayleigh  theory.  The  volume  seattering  funetion 
is  given  from  [3] 

A(«)  =  -  D’  (»’  +  2)"(1  +  cos«) 

lo/i 

where  v  is  the  thermal  eompressivity,  k  is  the  Boltzmann  eonstant,  n  is  the  refraetive 
index  and  T  is  the  absolute  temperature.  This  theory  deelares  that  seattering  also  depends 
on  Temperature  and  pressure  and  is  more  suitable  for  liquids  seattering. 


c.  Mie  Theory 

The  total  seattered  radiation  in  the  0  direetion  from  a  randomly  polarized 
beam  will  be  [3]. 
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I  =  —  [ (/j  +4)sin6'(i6’  /  =  —  [ {\  +i^)?,m6d6 
An ^  ^ 


An ' 


0  ''''  0 
where  \  is  the  intensity  scattered  in  the  0  direction  by  a  single  isotropic  sphere  from  an 
incident  polarized  beam  with  electric  vector  perpendicular  to  the  plane  of  observation  and 
4  is  the  intensity  in  the  plane  of  observation.  For  spherical  particles,  intensity  includes 

Riccati — Bessel  functions  and  size  parameters  ,  where  k^=2n  !  A  and  a  is  the  particle 
radius. 


4,  Volume  Scattering  Function 

Since  it  plays  a  key  role  in  the  theory  of  radiative  transfer,  the  volume  scattering 
function  is  the  most  important  inherent  property  for  the  understanding  of  optical 
oceanography.  Volume  scattering  P(0)  describes  the  probability  that  a  photon  will  be 
scattered  at  an  angle  0  (measured  from  the  forward  direction)  after  scattering  in  a  volume 
of  water.  The  larger  the  value  of  P(0)  is,  the  faster  the  beam  disperses. 


Figure  4.  Comparison  between  different  normalized  volume  scattering  functions  of 
ocean  water  as  a  function  of  the  scattering  angle  0.  (From  [3]) 
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5,  Dispersion  of  Scattering 

Scattering  has  a  wavelength  selectivity,  which  is  mainly  attributed  to  three 
factors:  Scattering  by  water  itself;  scattering  by  small-sized  particles  ;  and  absorption  by 
particles.  This  selectivity  can  be  seen  from  equation  (3.1)  and  the  factor  for  Rayleigh 
scattering. 


6,  Multiple  Scattering 

Multiple  scattering  is  the  irradiation  of  the  each  water  volume  element,  not  only 
by  the  original  beam,  but  from  light  scattered  from  previous  elements  as  well.  Thus,  the 
light  is  scattered  several  times.  This  effect  is  highly  dependent  on  the  particle 
concentration,  the  size  of  the  volume  and  the  attenuation  factor.  For  small  volumes  and 
low  particle  concentration,  multiple  scattering  is  considered  negligible. 

E,  ATTENUATION  OF  SEA  WATER 

The  combined  action  of  absorption  and  scattering  results  in  the  attenuation 
process  inside  the  sea  water.  The  important  parameter  describing  this  phenomenon  is  the 
total  attenuation  coefficient  c,  which  is  the  sum  of  the  scattering  and  the  absorption 
coefficients.  Attenuation  is  due  to  particles  and  due  to  absorption  of  dissolved  organic 
substances.  Sea  salt  has  little  effect  on  water  attenuation. 

Attenuation  is  primarily  an  absorption  that  depends  on  wave  length.  In  the  infra 
red  region,  scattering  can  be  completely  neglected  compared  to  absorption.  Generally, 
absorption  tends  to  increase  with  increasing  wavelength,  as  shown  in  Figure  5.  This 
relationship  is  attributed  to  several  infrared  absorption  bands. 

The  absorption  coefficient,  a(^),  is  a  measure  of  the  conversion  of  radiant  energy 
to  heat  and  chemical  energy.  It  is  numerically  equal  to  the  fraction  of  energy  absorbed 
from  a  light  beam  per  unit  of  distance  traveled  in  an  absorbing  medium.  Its  dependence 
on  the  wavelength  is  shown  in  Figure  5. 
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Figure  5.  Absorption  coefficient  a  as  function  of  wavelength  (From  [1]) 

When  dealing  with  attenuation  from  particles,  however,  the  wave  length 
dependence  is  different.  Attenuation  tends  to  increase  exponentially  at  shorter 
wavelengths,  as  shown  in  Figure  6.  Even  if  scattering  demonstrates  wave  length 
selectivity  as  discussed  before,  the  main  contribution  to  the  selective  effect  is  due  to 
absorption. 


itO  tCO  600  VOO  —  eOO 


Figure  6.  Particle  attenuation  coefficient  as  a  function  of  the  wavelength  (From  [3]) 
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III.  SIMULATING  THE  OCEAN  SURFACE 


This  chapter  introduces  some  of  the  background  knowledge  and  concepts  required 
for  understanding  wind-generated  ocean  surface  models. 

A,  THEORETICAL  APPROACH 

A  statistical  description  of  wind-generated  surface  waves  is  very  important  for 
examining  the  effects  of  the  laser  propagation  through  the  ocean  surface.  Much  work  has 
been  devoted  to  the  study  of  the  statistics  of  the  ocean  surface.  In  1953,  Cox  and  Munk 
studied  sun  glitter  on  the  ocean  surface  using  aerial  photography.  Oceanographers 
Pierson  and  Moskowitz  developed  spectral  density  functions  and  analytical  solutions  for 
the  height  variance  of  the  ocean  surface  versus  wind  speed.  Their  model  was  later 
improved  by  the  North  Sea  Wave  Project  (JONSWAP)  spectrum.  Elfouhaily,  Katsaros 
and  Chapron  refined  this  model  to  enhance  wind  effects  for  small  scale  waves. 

Regardless  of  the  means  of  construction,  a  wave  spectrum  must  agree  with 
observations,  as  well  as  satisfy  certain  criteria.  These  criteria  include  a  dependence  on 
fetch  conditions.  In  the  high-frequency  region,  the  wind  -  dependent  mean-squared  slope 
should  agree  with  observations  made  by  Cox  and  Munk. 

1.  Cox  and  Munk  Results 

Since  these  measurements  are  the  basis  for  evaluating  various  ocean  spectrum 
models,  the  results  obtained  by  Cox  and  Munk  are  presented  first. 

Combining  aerial  photography  of  the  sun’s  glitter  on  the  sea  surface  in  Hawaiian 
sea  and  measurements  of  wind  from  a  vessel  provided  the  following  results  regarding  the 
slope  distribution  of  the  surface.  Cox  and  Munk  showed  that  the  slope  variance  (mean- 
squared  slope),  regardless  of  direction,  increases  linearly  with  wind  speed  [4]. 
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O'"  =0.003  +  1.92x10  5 

erf  =3.16x10 

erf +a-f  =0.003  +  5.12x10  '6/i55 

where  erf  and  erf  eorrespond  to  the  eross  and  up  wind  eomponents  of  the  mean-squared 

slope.  5  represents  the  wind  speed  19.5  m  above  the  oeean  surfaee.  The  above  linear 

dependenee  is  seen  in  Figure  7.  The  Cox  Munk  results  have  beeome  the  referenee  eriteria 
for  testing  the  validity  of  an  oeean  wave  slope  measurement. 


Figure  7.  Mean-squared  slope  and  their  eomponents  as  a  funetion  of  the  wind  speed. 

(From  [4]) 
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2,  Pierson-Moscowitz  Spectrum 

The  Pierson-Moscowitz  Spectrum  has  the  advantage  of  being  simple,  reasonably 
close  to  actual  ocean  spectrum,  particularly  for  long  waves.  It  is  a  good  starting  point  for 
examining  the  various  models.  It  assumes  a  fully  developed  sea  (a  sea  produced  by  winds 
that  blow  steadily  for  a  long  period  of  time  over  a  large  area),  meaning  that  a  large  (>18 
km)  fetch  exists.  One  critical  benefit  is  that  the  Pierson-Moscowitz  Spectrum  has  a 
known  analytical  variance  and  surface  slope  spectrum.  The  ocean  surface  wave  power 
spectral  density  for  a  fully  developed  ocean  is  given  by  [5] 

S(ffl)  =  ^exp[-/?(^r]  (2.1) 

(O  CO 


In  this  expression,  a  =  0.0081,  P  =  0.74,  (Oq=—(Oq=  —  ,  where  g  is  the 

acceleration  of  gravity  and  U  is  the  wind  speed  at  19.5  m  above  the  surface.  The  angular 
frequency  (Z»  =  2;r/ ,  where  f  is  the  wave  frequency  in  Hz,  in  Eq.  (1.1)  is  defined  for 
positive  frequencies  only. 

To  depict  the  ocean  surface  in  spatial  coordinates,  the  Pierson-Moscowitz 
Spectrum  has  to  be  rewritten  in  terms  of  the  wave  number  k.  Using  the  deep  water 
dispersion  relationship 

o>^=gk  (2.2) 

and 

S(k)  =  S(a)^  (2.3) 

dk 

the  wave  number  spectrum  is  then  in  one  dimension 

^(0  =  ^exp[-(|U)]  ,2.4) 
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Equation  (1.4)  is  defined  for  positive  k  only.  Integrating  either  Eq.  (2.1)  or  Eq. 
(2.4)  over  co  or  k  gives  the  varianee  of  the  oeean  surfaee  wave  as 

^2  ^  ^  0.00274^^  (2.5) 

In  two  dimensions,  a  direetivity  faetor  is  needed  to  eompensate  for  the  laek  of 
isotropy  with  respeet  to  waves  in  the  direetion  of  wind  and  perpendieular  to  the  wind 
direetion.  The  two  dimensional  direetivity  faetor  is  given  by 

D{k,(p)  =  ^^  (2.6) 

n 

The  2D  speetrum  has  the  form  S(k)  =  ^jexp[-(^^)]D(A:,^)  (2.7) 

The  elevation  speetrum  for  wind  speeds  2-20  is  presented  in  Figure  8. 


Pierson  -  Moskowitz  Spectrum  for  2-20m/s  windspeeds 


Figure  8.  Pierson-Moskowitz  1-D  Elevation  Speetrum  for  different  wind  speeds  as  a 

function  of  the  wave  number  k 
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3,  Simulating  Ocean  Waves  with  Pierson-Moskowitz  Spectrum 

To  simulate  ocean  surface  waves,  a  Fourier  Transform  of  Equation  (2.4)  will  give 
the  autocorrelation  function  of  the  ocean  surface  waves.  This  method  does  not  give  the 
amplitude  versus  spatial  coordinates.  One  way  to  perform  the  numerical  simulation  is  to 
filter  Gaussian  random  noise  in  spatial  coordinates  using  the  square  root  of  Equation  (2.4) 
and  inverse  Eourier  Transforming  this  filtered  spectrum  to  spatial  coordinates.  A  more 
detailed  synopsis  of  the  procedure  required  can  be  presented  by  the  following  steps; 

•  Start  with  an  1-D  array  of  Random  Gaussian 

•  Fourier  Transform  this  Array  in  Frequency  Coordinates 

•  Filter  this  Spectrum  with  the  Pierson-Moscowitz  Spectrum 

•  Inverse  Fourier  Transform  Back  to  X  Coordinates 

•  Scale  the  Results  to  Correct  for  1/N  Factors 

The  Matlab  code  for  generating  a  numerical  simulation  of  the  ocean  surface  using 
the  Pierson-Moskowitz  Spectrum  is  code  A  in  Appendix  B. 


Noise  Filtered  Pierson-Moskowitz  Ocean  Spectrum 


Eigure  9.  Noise  filtered,  ID,  Pierson-Moskowitz  Spectrum  as  a  function  of  distance 
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4,  PM  Spectrum  Mean-Squared  Slope 

One  problem  with  the  Pierson-Moskowitz  Speetrum  is  that  it  underestimates  the 
high  frequeney  components  of  the  spectrum.  These  short  waves,  also  known  as  capillary 
waves,  have  an  important  contribution  to  the  mean-squared  slope  of  the  surface.  This 
behavior  can  be  observed  from  Figure  10,  where  the  dependence  is  obvious. 

As  described  previously,  a  wave  spectrum  must  be  capable  of  providing 
agreement  with  the  Cox  and  Munk  results  regarding  the  mean-squared  slope.  The 
Pierson-Moskowitz  Spectrum  diverges  logarithmically  as  the  wave  number  gets  bigger. 
This  is  due  to  the  dependence.  Analytically,  this  can  also  be  confirmed.  The  mean- 
squared  slope  is  obtained  by  integrating  the  spectrum  Eq.  (2.4)  over  all  wave  numbers  k 
as  follows  [6] 

si  =  =  (4.1) 

where  is  the  cut  off  wave  number. 

As  shown  in  Equation  (2.4),  S(k)  has  an  -independence.  Substituting  in  Eq.(4.1) 

K 

1  f  1 

inside  the  integral  retains  the  —  factor.  Our  initial  equation  has  now  become  . 

K 

This  integral  finally  gives  us  the  mean-squared  slope  In  . 

As  wave  number  k  increases,  the  mean-squared  slope  will  diverge,  reaching 
infinite  values.  This  is  not  physically  observed.  Ocean  waves  cannot  have  infinite  slope. 
Using  the  numerical  simulation  of  the  Pierson-Moskowitz  Spectrum,  the  results  are  the 
same.  Increasing  the  number  of  data  points  in  the  simulation,  the  mean-squared  slope 
constantly  increases.  This  is  without  saturating  to  a  specific  value.  Although  this 
spectrum  is  suitable  for  examining  the  low  frequency  region  (large  waves)  and  has  an 
analytical  form  to  calculate  the  mean-squared  slope,  it  does  not  fall  off  fast  enough  at 
high  frequencies. 


18 


5. 


JONSWAP  Model 


The  JONSWAP  spectrum  was  developed  for  the  North  Sea  oil  fields.  It  assumes 
that  a  wave  spectrum  is  never  fully  developed,  but  keeps  developing  through  non-linear, 
wave  interactions  for  long  distances  and  for  long  periods  of  time.  To  apply  this  theory, 
the  JONSWAP  spectrum  has  a  spectral  peak  enhancement  factor  Jp  included  in  the 

Pierson-Moskowitz  Spectrum.  This  factor  improves  the  spectrum’s  accuracy  at  low 
frequencies  by  including  the  wave  age.  The  problem,  however,  with  the  high  frequency 
waves  remains.  The  JONSWAP  factor  is  given  by  [7] 

Jp=Y^,  (6.1) 

where 

Y=  1.7  for  0.84  <  O  <  1  and 
Y  =  1 .7  +  6  In  (O)  for  1  <  n  <  5 


and 


r  =  exp[ 


2o-" 


(6.2) 


where  a  =  0.08  ( 1  +  4Q  ^ ).  The  inverse  wave  age  Q  is  defined  as 

Q  =  0.84tanh^-'' 


/  \0.4 

'  X  ^ 


(6.3) 


Aq  is  a  fetch  constant  and  has  the  value  =22000  m  and  the  non  -dimensional  fetch 


X  =  ,  where  x  is  the  dimensional  fetch,  x  =  10^  m. 


6/, 


10 
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An  improved  spectrum  that  includes  more  data  at  the  high  frequency  region  is, 
consequently,  more  suitable  for  this  study’s  purpose.  This  improved  spectrum  will  be 
presented  next. 


6,  Elfouhaily  Model 

The  model  proposed  by  Elfouahily  [8]  is  valid  for  a  larger  range  of  frequencies 
and  addresses  the  higher  frequency  waves.  Furthermore,  it  agrees  both  with  Cox  Munk 
data  and  the  Pierson-Moskowitz  analytical  results.  The  Elfouahily  model  is  more 
satisfactory  for  this  study’s  approach.  In  1-D,  the  Elfouahily  model  has  two  spectral 
regions  as  follows  [8] 

S(»r)  =  T[s+ii]  (7,1) 

K 

where  and  S,  express  the  high  and  the  low  frequency  curvature  spectrum.  The  total 

curvature  spectrum  is  expressed  as  S  =  /c^S  in  order  to  remove  the  dependence  of 
the  elevation  spectrum,  observed  in  the  previously  discussed  models.  The  rate  of  change 
of  the  ocean  surface  slope  can  be  approached  by  plotting  the  spectrum  as  a  function  of  the 
wave  number  k  as  presented  in  Figure  10.  The  effect  of  high  frequency  waves  is 
included.  After  the  value  k=300  rad/m,  the  spectrum  falls  off  very  fast.  This  prevents  the 
integral  of  the  spectrum  from  diverging. 
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Figure  10.  The  Elfouhaily  ID  Elevation  speetrum  S(k)  and  the  eorresponding 
curvature  spectrum  B(k)  as  a  function  of  the  wavenumber  k  for  different  wind 

speeds.  (After  [8]) 


The  low  frequency  term  (long  wave  curvature  spectrum)  is  defined  as 


Q 


exp[ — j=( - 1)] 


2  2 


Vio 


(7.2) 


where  ci^  is  the  Phillips-Kitaigorodskii  equilibrium  range  parameter  for  long  waves  and  is 
given  by  =  O.OOSyQ,  with  D  the  dimensionless  inverse  wave  age  parameter  defined  in 
Equation  (6.3).  The  phase  speed  at  spectral  peak  is  defined  as 


c 


p 


(7.3) 


where  the  wave  peak  wave  number  k  is  given  by  k  =Kq-  O? 


2 

and  kr,  =  — ^ 
0  ^2 


Lp^  is  the  Pierson-Moskowitz  spectral  shape  and  is  given  [8]  by 
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(7.4) 


The  high  frequency  term  in  Eq.  (6.1)  is  given  by  [8] 

B,=^-^-exp[-i(— -1)^]  (7.5) 

2  c  4  /r„, 

m 

where  is  the  Capillary  peak  celerity  and  has  the  constant  value  c^=0.  23  m/sec.  is 
the  Capillary  peak  wave  number  and  has  the  value  =36 1.4/m. 


To  define  the  parameter  (analogous  to  equilibrium  range  parameter  in  the  low 
frequency  spectrum),  the  friction  velocity  u  has  to  be  defined  as  [8] 


u  =  ■ 


0.42 -f/, 


10 


In 


10 


(7.6) 


where  Zg  is  the  roughness  length  and  is  given  by 


z,  =3.7-10*’ .(til/). 


0.9 


(7.7) 


The  above  equation  found  in  [8]  cannot  be  correct.  Checking  the  units  for  the 
roughness  length,  it  is  concluded  that  the  correct  form  has  to  be 


U  TJ 

5  /‘-^lO  A  /‘-^lO'.O.O 


Zg  =3.7-10^  •(^^)-(^)' 
g  c 


(7.8) 


in  order  to  get  meters  as  a  result.  Now  the  parameter  is  defined  as 


* 

u  * 

a^=  0.01-[l  +  ln( — )]  for  u  <c^ 


a„  =  0.01[l  +  31n(— )  for  u  >c„ 

As  seen  by  comparing  equations  (7.2)  and  (7.5),  the  Pierson-Moskowitz  spectral 

shape  is  only  contained  in  the  low  frequency  curvature  spectrum.  To  address  the  high 
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frequency  waves,  however,  this  low  frequency  roll-off  term  has  to  be  applied  in  both 
spectral  terms.  Without  this  term,  the  Elfouhaily  model  would  diverge  at  low  frequencies. 
This  is  because  the  high  frequency  curvature  spectrum  interacts  with  all  wave 

numbers.  Without  including  a  Pierson-Moskowitz  roll-off  the  spectrum  cannot  go  to  zero. 
Since  the  waves  would  have  infinite  height  at  low  frequencies,  this  cannot  be  a  physical 
consequence. 

B,  COMPUTER  SIMULATION 

After  deciding  upon  a  spectrum  model,  this  study’s  researchers  determined  that 
the  mean-squared  slope  of  the  Elfouhaily  Spectrum  must  be  compared  with  the  Cox 
Munk  results.  This  model  will  be  presented  in  the  third  part  of  this  thesis,  where  the 
Snell’s  Eaw  will  be  applied  in  an  ocean  surface. 

1,  Agreement  with  Observations 

Although  the  Elfouhaily  Spectrum  provides  a  result  with  a  physical  meaning,  it 
was  found  that  it  did  not  agree  with  the  Cox-Munk  results,  as  indicated  in  Eigure  11. 
Eurthermore,  it  was  observed  that  the  two  methods  used  for  slope  calculations  in  the  code 
B  seemed  to  deviate  from  each  other,  giving  slightly  different  results.  The  third 
observation,  which  seems  not  to  be  in  order,  is  a  discontinuity  that  exists  in  the  graph 
linearity,  especially  in  high  winds.  Eigure  1 1  shows  the  mean-squared  slope  versus  wind 
speed,  both  for  slope  1  and  slope  2. 
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slope  vs  wind  speed  for  N=  2  15 


wind  speed 

Figure  1 1 .  Slope  1  and  2  of  the  Elfouhaily  Spectrum  as  a  function  of  the  wind  speed 

compared  to  the  Cox  Munk  results. 

2,  Dependence  on  the  High  Frequency  Cut  Off 

Since  the  results  of  this  study’s  model  were  not  satisfactory,  changes  were 
applied.  A  first  look  indicates  that  the  model  works  well  in  low  wind  speeds,  but,  as  the 
wind  increases,  there  is  a  significant  deviation  from  Cox  and  Munk  results.  To  make  the 
simulation  more  realistic,  the  high  frequency  limit  was  first  increased  by  increasing  the 
number  of  data  points.  The  difference  between  the  two  slopes  disappeared.  In  Figurel2, 
this  result  is  presented  using  2^°  data  points  for  the  FFT  spectrum. 
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slope  vs  wind  speed  for  N=  2  20 


Figure  12.  Slopes  1  and  2  of  the  Elfouahily  Spectrum  as  a  function  of  the  wind  speed 

compared  to  the  Cox  Munk  results 


3,  Changing  the  Capillary  Peak  Wave  Number 

The  capillary  peak  wave  number  is  a  function  of  the  water  density  p,  the 

gravitational  constant  g  and  the  ocean  surface  tension  T.  It  determines  the  sharp  roll  off 
of  the  Elfouahily  Spectrum  at  high  frequencies.  It  is  given  by 


Pg 

T 


and,  as  stated  before,  has  the  value  of  361 .4(  rad  /m). 

Eield  measurements  of  intermediate  scale  waves  by  Flwang  [9],  however,  showed 
different  results.  The  observed  peak  occurs  at  much  lower  wave  numbers  than  the  above 
constant.  Eurthermore,  is  not  constant,  but  depends  on  wind  speed.  This  data  is 
presented  in  Eigure  14. 

The  area  under  the  curve  in  a  range  of  wave  numbers  corresponds  to  the  mean- 
squared  slope,  since  = /r^5'(x') .  The  calculated  slopes  are  in  agreement  with  those 
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obtained  from  [4],  Observing  the  peak  of  the  distribution  of  the  mean-squared  slopes 
represented  by  B{k)  allows  the  estimation  of  a  more  realistic  value  for  capillary  peak 
wave  number  in  the  simulation  code. 


Figure  13.  Curvature  spectrum  measured  by  Hwang  as  a  function  of  the  ocean  wave 

number  k  (From  [9]) 

As  can  be  seen,  the  peak  of  the  spectrum  for  various  wind  speeds  appears,  at  a 
maximum,  in  ranges  from  20  rad/m  to  200  rad/m.  Making  different  trials  in  the 
simulation  code  and  comparing  the  graphic  results  with  the  corresponding  Cox  Munk 
data,  these  researchers  selected  the  value  of  160  rad/m  for  Km.  The  mean-squared  slope 
graph  can  be  seen  in  Figure  14.  It  can  now  be  seen  that,  for  smaller  wind  speeds,  the 
corresponding  data  almost  overlaps,  demonstrating  that  the  Elfouhaily  model  spectrum  is 
appropriate  for  these  wind  speeds.  For  stronger  winds,  however,  there  seems  to  be  only  a 
small  deviation. 
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slope  vs  wind  speed  for  N=  2  20 


wind  speed 

Figure  14.  Slopes  1  and  2  of  the  Elfoyhaily  Spectrum  compared  to  the  Cox  Munk 

results  for  Km  =160 

4,  Capillary  Peak  Wave  Number  as  a  Function  of  Wind  Speed 

Taking  again  into  consideration  Hwang’s  results,  these  researchers  considered  the 
relation  between  capillary  peak  wave  number  and  speed.  Hwang  showed  that  shifts 

upward  with  increasing  wind  speed.  Analyzing  the  data  that  corresponded  to  the  peak  of 
the  spectrum  for  each  used  wind  speed  from  figure,  the  results  ended  in  a  linear 
relationship  between  the  two  parameters 


=0.077872^-31.797  (4.1) 


Substituting  this  relationship  in  this  study’s  code  and  running  again  the 
simulation,  the  results  are  shown  in  Figure  15.  A  good  agreement  with  the  desired  results 
appears  in  the  higher  wind  speeds  and  an  approach  that  is  satisfactory  at  lower  speeds. 
Furthermore,  there  is  no  discontinuity  in  the  graph’s  linearity,  something  that  was 
common  for  a  fixed  capillary  wave  spectral  peak. 
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The  poor  agreement,  seen  in  Figure  15  at  low  frequencies,  suggests  that  a 
quadratic  dependence  on  U  may  be  more  appropriate  than  the  linear  dependence  used  in 
equation  (4.1). 
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Figure  15.  Slopes  1  and  2  of  the  Elfoyhaily  Spectrum  as  a  function  of  the  wind  speed 
compared  to  the  Cox  Munk  results,  using  the  relationship  Km  =  12.69  U  -31.797 

C.  CHAPTER  SUMMARY 

In  this  chapter,  two  different  ocean  spectrums  were  used  to  generate  a 
mathematical  model  for  a  wind-generated  ocean  surface.  Initially,  the  Pierson-Moskowitz 
Spectrum  was  examined,  but  rejected  since  the  corresponding  mean-squared  slope 
diverged  at  high  frequencies.  Next,  the  Elfouhaily  Spectrum  was  tested.  Although  more 
complicated,  it  did  agree  with  physical  results,  having  a  non  diverging  mean-squared 
slope.  The  agreement  with  measured  observations  was  then  evaluated.  After  altering 
specific  parameters  in  the  mathematical  approach  and  computer  simulation,  an  ocean 
surface  model  was  constructed  that  agreed  with  Cox  Munk  observations,  a  fact  that 
reinforces  its  validity  and  suitability  for  this  study’s  purpose. 
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IV.  SIMULATION  AT  THE  AIR  WATER  INTERFACE 


A,  GEOMETRICAL  OPTICS’  APPROACH 

This  section  develops  a  geometrical  optics’  model  to  examine  the  effects  of  light 
entering  or  leaving  a  water  surface.  The  ocean  model  discussed  in  Chapter  III  was  used. 
To  generate  a  physically  accurate  ocean  surface  from  this  model,  this  study  followed  the 
same  procedure  as  described  in  paragraph  II-3.  White  noise  was  fdtered  with  the  wave 
number  spectrum.  The  final  outcome  was  the  ocean  slope  distribution.  The  incoming 
light  ray  interacts  with  this  slope,  refracting  according  to  the  Snell’s  Law  and  exiting  with 
a  new  angle.  The  corresponding  code  is  in  the  Appendix. 

1.  Snell’s  Law 

The  simulation  geometry  is  as  follows:  A  light  ray/j  is  incident  on  the  water 
surface  with  corresponding  angle  /j  to  the  normal.  The  ray  exits  with  angle  •  The  local 
ocean  slope  is  R.  By  Snell’s  Law  we  have 


n,  sin  y-j  =  sin 


(1.1) 


where  n,  =  1.3  3,^2  =1  are  the  corresponding  indices  of  refraction  for  water  and  air. 
Figure  16  represents  the  geometry  of  the  simulation. 
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Figure  16.  Geometry  Simulation  for  applying  Snell’s  Law  in  an  ocean  surface  with 

slope  distribution  R  (After  [10]) 


Solving  for  ,  we  get  [1 1], 

Tl 

=  sin  '  [—  sin(/2  +  R)]  -  R 


(1.2) 


This  relationship  applies  to  the  laser  light  distribution  just  beneath  the  boundary 
and  is  valid  for  light  both  entering  and  exiting  the  ocean  surface.  The  entering  and  exiting 
angles  are  always  defined  with  respect  to  the  vertical. 

2.  Assumptions 

The  scattering  and  absorption  effect  will  be  initially  neglected  and  examined  in 
the  following  chapter.  Scattering  includes  both  the  initial  interference  and  the  following 
multi-scattering  phenomena.  Wave  obscuration  of  the  optical  rays  is  also  neglected. 
Propagation  is  assumed  to  be  vertical  and  the  ocean  angular  variance  R  is  small. 
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3,  Analytical  Solutions 

From  [11],  the  variance  of  the  ray  angle  after  penetrating  the  sea  surface  is  given 
approximately  by 


var[yj]=^ 


1 


var[i?] 


(3.1) 


and 


var[i?] 


(3.2). 


Calculating  the 


ydRj 


from  equation  [1.2]  with  respect  to  R  results  in 


dy, 

dR 


cos(/  +  (R)) 

1 

2 

tl 

h- 

sin"(/2+(^)) 

V 

) 

-1 


(3.3) 


Squaring  the  result  and  multiplying  by  the  variance  of  the  local  slope  R  gives  a  better 
approximation  of  the  variance  of  the  refracted  angles. 


var 


ydR  y 


var[R] 


(3.4) 


4.  Simulation 

The  simulation  was  written  in  Matlab  and  code  C  and  it  appears  in  Appendix. 
Local  slope  R  is  determined  from  the  Elfouhaily  Spectrum  and  was  calculated  with 
simulation  code.  The  output  is  shown  in  Figure  17,  which  demonstrates  the  ratio  of 
variances  versus  the  incident  angle. 
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Figure  17.  Variance  of  ratio  of  exiting  angle  and  local  slope  versus  incident  angle 

B.  CHAPTER  SUMMARY 

This  chapter  used  Snell’s  Law  to  determine  the  distribution  of  refracted  waves 
versus  wind  speed  to  the  ocean  surface  model  generated  in  the  previous  chapter.  A 
geometric  optics’  approach  was  used.  Scattering  and  absorption  during  light’s  penetration 
in  the  water  body  were  ignored.  Examining  the  ratio  of  variances  of  the  exiting  angle  and 
the  local  slope,  it  was  shown  that  the  refractive  effects  of  the  ocean  surface  were 
relatively  small. 
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V.  THEORY  OF  RADIATIVE  TRANSFER  IN  THE  SEA 


The  complex  behavior  of  the  photons  inside  a  water  body,  from  the  simultaneous 
effects  of  absorption  and  scattering  and  the  shape  of  the  volume  scattering  function, 
prevents  the  establishment  of  an  analytical  relationship  between  the  apparent  and  the 
inherent  optical  properties.  Monte  Carlo  Computer  simulation  techniques  can  provide  a 
description  of  the  scattered  underwater  light  field,  given  the  measured  optical  properties 
of  the  water  and  the  characteristics  of  the  incident  radiant  flux. 

A,  MONTE  CARLO  SIMULATION  OF  LIGHT  PENETRATION  INTO 

NATURAL  WATERS 

Although  the  sea’s  inherent  optical  properties  that  determine  the  light  behavior 
inside  water  may  be  known,  there  is  no  analytical  means  to  predict  the  photons’ 
movement  in  a  scattering  medium.  Thus,  the  underwater  light  distribution  field  cannot  be 
determined  directly. 

These  optical  properties,  however,  do  provide  useful  information.  They  inform  us 
of  all  the  possible  ways  a  photon  can  move,  the  possible  interactions  with  the  scattering 
medium  and  the  corresponding  probabilities  of  the  above.  A  photon  can  be  absorbed  in  a 
specific  path  length  or  scattered  in  a  specific  path  length  and  at  a  specific  angle.  Using 
suitable  computer  simulation  techniques,  for  given  values  of  inherent  properties,  the 
corresponding  apparent  ones  and  other  parameters  of  the  underwater  light  field  have  been 
calculated  by  Kirk. 

This  is  the  idea  behind  the  Monte  Carlo  simulation  where  the  movement  of  a  large 
number  of  photons  is  calculated  in  the  computer  and  the  average  photon  flux  is 
calculated.  The  code  implemented  by  Kirk  using  Monte  Carlo  simulation  was  used  for 
the  scope  of  this  thesis.  A  brief  synopsis  is  provided  next. 
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1. 


Kirk  Code 


Reflection  from  the  surface  was  ignored.  The  water  surface  was  considered  flat 
and  only  the  direct  light  beam  was  taken  into  consideration.  The  water  body  was  assumed 
to  have  finite  depth.  The  objective  was  to  estimate  the  nature  of  light  field  at  each  depth. 

Just  below  the  surface  and  after  refraction,  a  photon’s  trajectory  is  determined  by 
its  angle  a  with  respect  to  the  horizontal  given  by 

«  =  arccos[(cos  «*)/!. 33]  ,  (1.1) 

where  a  is  the  light  altitude.  The  probability  of  the  photon  moving  at  a  specific  distance 
before  its  first  interaction  is  determined  by  the  attenuation  coefficient  c.  A  random 
number  r  between  0  and  1  is  selected  and  the  corresponding  path  length  is  [12] 

P  =  -(-)ln(l-r)  (1.2) 

c 

Psin(  a )  gives  the  corresponding  depth  .  A  new  random  number  decides  the  nature  of  the 
interaction.  For  absorption,  a  new  photon  enters  the  water  and  the  procedure  continues. 
Scattering  requires  a  new  scattering  angle  0  (angle  between  photon’s  trajectory  before 
and  after  scattering)  and  rotation  angle  cp  (angle  of  rotation  of  the  new  trajectory  around 
the  previous  one).  The  probability  that  a  photon  may  be  scattered  at  an  angle  0  is  given 
from  the  volume  scattering  phase  function  P(0).  A  new  random  number  is  assigned  a 
value  of  scattering  angle  0  from  the  cumulative  distribution.  The  water  had  a  normalized 
volume  scattering  function  P(0),  which  is  the  same  as  the  experimental  measurements  of 
Petzold  in  San  Diego  Harbor.  A  subsequent  random  number  provides  the  rotation  angle 
(p.  The  new  angle  a'of  the  photon  relative  to  the  horizontal  is  now  [12] 

sina'  =  sinacos6>  +  cosacos^  (1.3) 

where  a  is  the  angle  calculated  in  the  previous  step.  Before  the  next  interaction  and  the 
corresponding  depth,  a  new  random  number  determines  the  path  length. 

When  the  movements  of  all  photons  are  calculated,  the  total  photon  flux  at  each 
depth  provides  a  measure  of  the  downward  or  upward  radiance.  The  other  quantities 
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provided  as  outcome  from  this  program  are  the  average  vertical  attenuation  coefficient, 
average  cosine  and  radiance  distribution  in  each  depth  ,  the  x,  y  coordinates  and  the  0,  (p 
angles  for  each  photon. 

The  above  calculation  applies  to  monochromatic  radiation  only. 

2,  Testing  the  Kirk  Code 

To  validate  the  Monte  Carlo  code  and  to  provide  useful  and  reliable  results,  the 
results  should  agree  with  known  measurements,  as  in  Chapter  III,  which  compared 
Elfoyhaily  model  to  the  Cox  Munk  results.  The  code  results  were  initially  compared  to 
measured  values  of  irradiance  attenuation  coefficient  for  given  values  of  solar  elevation, 
wavelength,  scattering  and  absorption  coefficients.  The  input  values  used  were  the  same 
as  in  Table  10  from  [3].  Table  15  from  [3]  provided  the  scattering  coefficient  data.  The 
depth  interval  between  layers  of  Im  was  used  in  the  program  according  to  the 
measurements.  The  code  results  are  very  close  to  the  measured  ones,  as  indicated  in 
Tables  1  and  2. 

The  second  comparison  included  the  irradiance  attenuation  coefficient  from  [13]. 
This  was  for  the  same  given  scattering  parameters  and  scattering  function.  The  only 
difference  was  that  a  constant  solar  angle  of  90  degrees  and  a  depth  interval  of  0.1  m 
were  used.  Because  of  the  high  scattering  attenuation,  the  code  results  were  very 
accurate. 

Both  comparisons  proved  to  be  very  close  to  the  observed  values  and  so  validated 
the  code. 
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Wavelength 

(nm) 

Absorption 
coefficient  a 

(m-') 

Scattering 
coefficient  p 
(m  ‘) 

Irradianee 

attenuation 

coefficient 

Irradianee 

attenuation 

coefficient 

Ek 

(m') 

calculated  by 

Kirk  Code 

200 

3.07 

0.151 

3.14 

3.08 

300 

0.141 

0.0262 

0.154 

0.141 

400 

0.0171 

0.0076 

0.0209 

0.017 

500 

0.0257 

0.0029 

0.0271 

0.026 

600 

0.244 

0.0014 

0.245 

0.244 

700 

0.650 

0.0007 

0.650 

0.649 

800 

2.07 

0.0004 

2.07 

2.06 

Table  1 .  Comparison  between  observed  values  of  the  Irradianee  attenuation 


Region 

Solar 

elevation 

in 

degrees 

Wavelength 
k  (nm) 

Absorption 
coefficient 
a  (m“' ) 

Scattering 

coefficient 

P(m-') 

Irradianee 

coefficient 

a 

(m"') 

Irradianee 
coefficient 
calculated 
by  Kirk 
Code 

Mediterranean 

Sardenia 

72 

372 

0.058 

0.05 

0.07 

0.06 

Mediterranean 

Sardenia 

74 

633 

0.29 

0.05 

0.31 

0.30 

Baltie 

52 

372 

0.97 

0.2 

1.22 

1.11 
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Baltie 

48 

633 

0.33 

0.2 

0.42 

0.40 

Baltie 

52 

535 

0.09 

0.2 

0.12 

0.118 

Table  2.  Comparison  between  observed  values  of  the  Irradianee  attenuation 

eoeffleient  and  results  derived  from  the  Monte  Carlo  simulation 

B,  SIMULATION  RESULTS 


Running  the  Kirk  Code  for  various  input  values  of  both  absorption  and  seattering 
ooeffieients  provided  a  measure  of  the  effeet  of  these  parameters  for  light  distribution 
within  the  oeean.  The  oeean  water  body  was  divided  into  21  depth  layers.  After  the  21®* 
layer,  the  seattered  photon  distribution,  starting  from  the  initial  point  foeused  laser  beam, 
interaeted  with  the  loeal  slope  oeean  surfaee  from  the  Elfouhaily  model  using  Snell’s 
Law  .This  gave  the  final  distribution  of  the  photons  as  they  exit  the  surfaee. 

An  initial  photon  elevation  angle  of  90  degrees  and  a  depth  interval  of  10m  were 
used  throughout  the  simulation. 

The  original  eode  was  written  [12]  in  FORTRAN  77  and  was  translated  into 
FORTRAN  95  by  Professor  Walters.  This  translation  was  done  to  avoid  the  many  ‘go  to 
statements  and  baek  referenees.’ 

Three  eases  were  examined  independently.  In  the  first  ease,  named  I,  very  small 
values  of  the  absorption  and  seattering  eoeffieient  were  ehosen.  These  eorresponded  to 
negligible  effeets  of  the  absorption  and  seattering  phenomena.  Next,  in  ease  II,  higher 
values  were  seleeted  to  examine  the  relationship  of  the  above  parameters  with  the  light 
distribution  inside  the  water  when  the  above  phenomena  are  signifieant.  Finally,  the  last 
ease,  named  III,  eonsidered  intermediate  values  of  the  absorption  and  seattering 
eoeffieients.  This  was  mainly  to  test  and  eonfirm  that  there  is  a  logieal  and  expeeted 
relationship  between  the  light  distribution  and  the  laser  footprint  and  the  variation  of  the 
inherent  optieal  properties. 
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In  each  case,  two  3D  graphs  of  the  eleventh  and  twenty-first  layers  were 
presented.  The  last  one  is  the  most  important.  This  is  because  it  describes  the  light 
distribution  just  before  the  interaction  with  the  ocean  surface.  For  better  understanding  of 
the  laser  footprint  size  and  shape,  a  2D  graph  of  the  twenty-first  layer  is  also  provided. 

1.  Case  I  (Small  Extinction) 

For  this  case,  the  following  values  will  be  assigned  to  the  input  parameters; 

Absorption  coefficient  a  =  0.01  1/m 

Scattering  coefficient  b  =  0.01  1/m 

Depth  interval  =  1  m  for  a  total  depth  of  21  m 

Photon  angle  in  the  water  =  90  degrees 
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THETA  [rad] 


Figure  18.  Light  distribution  at  the  eleventh  layer  for  ease  1  for  10  m 


Figure  19.  Light  distribution  at  the  twenty  first  layer  for  ease  1  for  20  m 
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Figure  20.  2D  Light  distribution  at  the  twenty  first  layers  for  case  I  for  20  m 

2,  Case  II  (Higher  Extinction) 

For  this  case,  the  following  values  will  be  assigned  to  the  input  parameters; 

Absorption  coefficient  a  =  0.5  1/m 

Scattering  coefficient  b  =  0.5  1/m 

Depth  interval  =  0.5  m  for  a  total  depth  of  10.5  m 

Photon  angle  in  the  water  =  90  degrees 


2D  graph  of  the  laser  beam  footprint  at  the  21  st  layer 
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THETA  [rad] 


3D  graph  of  the  laser  beam  at  the  1 1  th  layer 
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Figure  21 .  Light  distribution  at  the  eleventh  layer  for  case  11  for  5  m 


3D  graph  of  the  laser  beam  at  the  21  st  layer 


Figure  22.  Light  distribution  at  the  twenty  first  layer  for  case  11  for  10  m 
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Figure  23.  2D  Light  distribution  at  the  twenty  first  layer  for  case  II  for  10  m 

3,  Case  III  (Moderate  Extinction) 

For  this  case,  the  following  values  will  be  assigned  to  the  input  parameters; 

Absorption  coefficient  a  =  0.1  1/m 

Scattering  coefficient  b  =  0. 1  1/m 

Depth  interval  =  1  m  for  a  total  depth  of  2 1  m 

Photon  angle  in  the  water  =  90  degrees 
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THETA  [rad] 


3D  graph  of  the  laser  beam  at  th  1 1  th  layer 
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Figure  24.  Light  distribution  at  the  eleventh  layer  for  ease  III  for  10  m 
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Figure  25.  Light  distribution  at  the  twenty  first  layer  for  ease  III  for  20  m 


3D  graph  of  the  laser  beam  at  the  21  st  layer 
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Figure  26.  2D  Light  distribution  at  the  twenty  first  layer  for  case  III  for  20  m 

C.  INTERACTION  WITH  THE  SEA  SURFACE 

This  chapter  evaluates  the  effect  of  the  sea  surface  on  the  light  distribution. 
Compared  to  Chapter  III,  the  scattering  and  absorption  influences  of  the  water  body 
before  the  interaction  were  included. 

The  light  distribution  just  before  the  interaction  was  calculated  in  paragraph  B  of 
this  chapter  using  both  3D  and  2D  graphs  of  the  twenty-first  layer  for  each  case,  as  well 
as  the  corresponding  data  contained  in  the  Kirk  Code. 

To  demonstrate  the  spread  of  the  laser  beam  as  it  exits  the  sea  surface,  a 
histogram  of  the  exiting  angle  will  be  given  for  case  III  in  comparison  to  the  histogram  of 
the  scattering  angle  just  before  the  interaction.  For  a  wind  speed  of  10  m/sec  the 
negligible  difference  between  the  two  plots  indicates  the  small  effect  of  the  ocean  surface 
refraction. 


2D  gtph  of  the  laser  beam  footprint  at  the  21  st  layer 
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Histogramm  of  the  scattering  angle  at  the  21  st  layer 
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Figure  27.  Scattering  angle  for  absorption  and  scattering  at  the  twenty  first  layer  for 

case  III.  Standard  deviation  =  0.2842 


Figure  28.  Angle  spreading  for  absorption,  scattering  and  refraction  for  case  III 

•  Standard  Deviation  =  0.3251 

D.  STANDARD  DEVIATION 

Applying  the  Kirk  Code,  a  matrix  table  for  standard  deviation  is  next  presented.  It 
indicates  the  effect  of  scattering  for  the  water  depth  of  twenty  one  meters.  As  it  can  be 
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seen  from  Table  3,  for  constant  absorption  coefficient,  standard  deviation  increases  with 
increasing  scattering  coefficient,  as  expected.  Finally,  scattering  appears  to  be  more 
important  after  the  value  of  0.5  (1/m). 

Table  3.  Standard  Deviation  of  the  scattering  angle  for  the  depth  of  21  m 


a  [l/m]\  b  [1/m] 

0.01 

0.05 

0.1 

0.5 

0.01 

0.2826 

0.2617 

0.2434 

0.2017 

0.05 

0.2947 

0.2832 

0.2548 

0.2218 

0.1 

0.3317 

0.3016 

0.2852 

0.2531 

0.5 

0.5044 

0.4686 

0.4443 

0.3600 

E,  CHAPTER  SUMMARY 

This  chapter,  which  is  the  core  of  the  thesis,  examines  initially  the  spreading  of 
the  laser  beam  propagating  inside  a  water  body  from  scattering  and  absorption.  It 
considers  theoretical  approach  of  the  Monte  Carlo  and  the  Kirk  Code  and  the  simulation 
results.  Examining  the  results,  it  can  be  seen  that  the  dominant  phenomenon  is  absorption 
compared  to  scattering.  The  upward  light  distribution  was  refracted  at  the  sea  surface, 
resulting  in  the  beam  pattern  exiting  the  water.  Comparing  the  histograms  of  the 
scattering  and  exiting  angles  it  was  found  that  the  main  effect  is  that  of  the  scattering 
itself.  Further,  compared  to  scattering  and  absorption,  the  refraction  at  the  surface  does 
not  substantially  alter  the  beam  pattern  created  as  indicated  by  the  corresponding  standard 
deviations,  which  are  very  close. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A,  THESIS  CONCLUSIONS 

A  laser  beam  fired  upward  inside  a  water  body  experienees  seattering  and 
absorption  that  expands  the  beam  of  light.  The  absorption  and  seattering  eoeffieients  and 
the  volume  seattering  funetion  determine  the  angular  spread.  These  parameters  vary  with 
wave  length  and  the  oeean  water  itself.  The  light  distribution  is  then  refraeted  at  the 
oeean  surfaee  giving  the  final  footprint  of  the  beam,  propagating  in  the  air. 

Among  the  two  phenomena,  known  together  as  extinetion,  the  absorption  is  the 
dominant  one.  Furthermore,  no  eonelusions  ean  be  made  based  independently  on 
absorption  only  or  seattering  only  eoeffieients.  This  is  beeause  their  aetion  is 
simultaneous  and  eombined.  Generally  the  spread  of  the  beam  inereases  as  absorption 
inereases.  The  photons  do  not  have  the  ability  to  propagate  far,  even  for  higher  values  of 
the  seattering  eoeffieient. 

Between  attenuation  inside  the  oeean  and  seattering  at  the  surfaee,  the  first  is 
more  important  and  mainly  defines  the  shape  and  the  size  of  the  laser  beam  footprint. 

B,  RECOMMENDATIONS 

Now  that  the  laser  beam  footprint  ean  be  approximated,  the  attenuation  due  to 
propagation  in  the  air  and  elouds  should  be  eonsidered.  The  objeetive  is  to  obtain  the 
final  light  distribution  at  the  desired  loeation. 

Aeeordingly,  when  the  laser  beam  is  fired  from  the  air  into  the  water,  the  above 
results  need  to  be  eonsidered. 

Another  issue,  apart  from  the  beam  spreading,  is  the  desired  laser  irradianee 
delivered  on  the  target  and  the  eorresponding  effeetiveness. 
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APPENDIX  A 


A,  FUNDAMENTAL  QUANTITIES 

1,  Radiant  Flux  {F} 

The  time  rate  of  flow  of  radiant  energy. 

W 

— ^  ,  where  Q  is  the  radiant  energy.  Measured  in  W. 
m 

2,  Radiance  {L} 

The  radiant  flux  per  unit  solid  angle  per  unit  projected  area  of  surface. 

,  d^F  W 

L  = - .  asured  in  — ^ — . 

dA  ■  cos  s  ■  dco  m  sr 

3,  Irradiance  {E} 

Irradiance  at  a  point  of  a  surface  is  the  radiant  flux  incident  on  an  infinitesimal 
element  of  a  surface  containing  the  specific  point,  divided  by  the  element  area. 

dF  W 

E  =  —  .  Measured  in  — ^  . 
dA  m 
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APPENDIX  B:  CODES 


A,  MATLAB  CODE  FOR  GENERATING  THE  PIERSON-MOSKOWITZ 
NOISE  FILTERED  SPECTRUM  AND  CALCULATING  THE  MEAN- 
SQUARED  SLOPE 


%  ocean  surface  and  variance 
%29  JULY  09 


%constants 
a  =  0.0081; 
b  =  0.74; 

g  =  9.81;  %gravity  acceleration  in  m/sec^2 
U  =  10;  %wind  speed  in  m/sec 
N  =  2.^15;  %  number  of  points 

%  parameters 
L  =  b. *g. ^2 . /U. ^4; 

dk  =  0.02;  %  linear  stepsize  of  wavenumber 

k  =  1 :N; 

k  =  dk*k; 

k2  =  N/2; 

k3  =  k2+l; 

k4  =  l:k2; 

k5  =  l:k2-l; 


%  PM  spectrum  calculation 

S  =  (a. / (2 . *k. "3) ) . * (exp (- (L. /k. "2) ) 
figure  (1) 
loglog ( k, S ) 

A  =  randn ( 1 , N) ; 
gaussian  numbers 
B  =  f ft  (A)  ; 

C (1 :N/2) =B (1 :N/2) . *sqrt (S (1 :N/2) ) ; 

PM  spectrum 
C(k3)=0; 

C ( k3+k5 ) =con j (C ( k3-k5 ) ) ; 

d  =  ifft(C); 

coordinates 

dx  =  (2*pi) / (N*dk) 

coordinates 

X  =  dx  * ( 1 : N ) ; 

d  =  d* (sqrt (2*pi) ) /sqrt (2*dx) ; 
figure  (2) 
plot (x, d) 

xlabel (' distance  in  m') 
ylabel (' height  in  m') 

%variance  calculation 


);  %  PM  spectrum  in  m^2/red/m 

%  PM  spectrum  graph 

%create  an  1-D  array  of  random 

%  noise  spectrum 

%  filter  of  noise  spectrum  with 

%  inverse  fft  to  distance 
%translate  k  into  x  distance 
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variance=  var(d) 

Var  =  dk. *sum (S) 

%slope  calculation 
slopel  =  dk*sum (k. ^2 . *S) 
y  =  diff ( (d) /dx) ; 

slope2  =  var(y) 

slopes  =  (a/4)*  expint  ( (b*g^2) . / (k (N) . ^2*U^4) ) 


B,  MATLAB  CODE  FOR  GENERATING  THE 
FILTERED  SPECTRUM  AND  CALCULATING 
SLOPE 


%  06  August 
%J0NSWAP  1-D  Spectrum 


in  m/sec^ 
in  m/ sec 
in  1/m 

in  m/ sec 

in  m 

%  parameters 
dk  =  0.02; 

k  =  1:N;  %  wavenumber 

k  =  dk*k; 

k2  =  N/2; 

k3  =  k2+l; 

k4  =  l:k2; 

k5  =  l:k2-l; 


%  general  costants 


g  = 

9.81; 

%gravity  acceleration 

Cm  = 

0.23; 

%capillary  peak  celerity 

Km  = 

160; 

%capillary  peak  wavenumber 

Xo  = 

22000; 

%  fetch  constant 

U  = 

10; 

%  wind  speed 

N  = 

2  .  "20; 

%  Number  of  points 

XI  = 

< 

o 
( — 1 

%  dimensional  fetch 

%  low  frequency  curvature  spectrum  B1  calculation 
%  special  costants 
Ko  =  g/U.^2  ; 

X  =  (g*Xl)/U^2;  % 

omega  =  0.84  . * (tanh ( (X/Xo)  . ^0 . 4) )  . ^  (-0 . 75) ; 

Kp  =  Ko . *  omega. ^2;  % 

wave  peak 

Ap  =0.006*  sqrt (omega) ;  % 

parameter 

Cp  =  sqrt  (g/Kp) ;  % 

spectral  peak 

c  =  sqrt  (g . /k)  ;  % 

Lpm  =  exp ( (-5/4) . * ( (Kp. /k) . ^2) ) ;  % 

s  =  0.08*(l  +  4.^  omega^-3)  ; 

GAMA  =  exp (- ( (sqrt (k. /Kp) -1)  . ^2) / (2*s^2) )  ; 


ELFOUHAILY  NOISE 
THE  MEAN-SQUARED 


non  dimensional  fetch 

wavenumber  of  gravity 

equilibrium  range 

phase  speed  at 

phase  speed 
PM  spectral  shape 
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if  omega  <1 

gama  =  1.7; 

else  gama  =  1 . 7+ln (omega) ; 
end 

Jp  =  gama.^GAMA;  %  JONSWAP  peak 

enhancement  factor 

Z  =  exp (  (-  omega/sqrt (10) ) . * (sqrt (k. /Kp) -1) ) ; 

B1  =  (Ap/2 ) * (Cp . /c) . *Lpm. *  Jp.*Z;  %low  frequency 

curvature  spectrum 

figure  (1) 

loglog ( k, B1 ) 

xlabel (' wavenumber  k') 

ylabel ( ' low  frequency  curvature  spectrum  B1 ' ) 

%high  frequency  curvature  spectrum  Bh  calculation 

Zo  =  3 . 7* (10^-5) * (U^2/g) * (U/Cp) ^0 . 9;  %  roughness  length 

ustar  =  0 . 42*U/log (10/Zo) ;  %friction  velocity 

if  ustar<Cm 

Am  =  0 . 01* (1+log (ustar/Cm) ) ; 

else  Am=  0 . 01* (l+3*log (ustar/Cm) ) ; 

end 

Bh  =  (Am/ 2 ) * (Cm. /c) . *exp (-0.25* ( ( k . /Km) -1 ) . ^2 ) . *Lpm;  %high  frequency 
curvature  spectrum 

figure  (2) 

loglog ( k, Bh) 

xlabel (' wavenumber  k') 

ylabel ('high  frequency  curvature  spectrum  Bh ' ) 

SPEC  =  (Bh+Bl) ./k.^3;  %  1-D  spectrum 

%  plotting  the  spectrum 
figure ( 3 ) 
loglog (k, SPEC) 
xlabel (' wavenumber  k') 
ylabel (' elevation  spectrum  SPEC(k) ') 

%  filtering  the  spectrum 

A  =  randn ( 1 , N) ; 
gaussian  numbers 
B  =  f ft  (A)  ; 

C (1 :N/2) =B (1 :N/2) . *sqrt (SPEC (1 :N/2) ) ; 

PM  spectrum 
C(k3)=0; 

C (k3+k5) =conj (C (k3-k5) ) ; 
d  =  ifft(C); 
coordinates 
dx  =  (2*pi) / (N*dk) ; 
coordinates 
X  =  dx  * ( 1 : N ) ; 


%create  an  1-D  array  of  random 
%  noise  spectrum 

%  filter  of  noise  spectrum  with 

%  inverse  fft  to  distance 
%translate  k  into  x  distance 


53 


d  =  d* (sqrt (2*pi) ) /sqrt (2*dx) ; 

%  plotting  the  noise  filtered  spectrum 
Figure  (4) 
plot (x, d) 

xlabel (' distance  in  m') 
ylabel (' height  in  m') 

title ('noise  filtered  SPEC  spectrum') 

%  variance  calculation 
variance  =  var(d) 

VAR  =  dk. *sum(SPEC) 

%  slope  calculation 
slopel  =  dk*sum (k. ^2 . *SPEC) 
y  =  diff ( (d) /dx) ; 

slope2  =  var(y) 

%  slopes  =  (a/4)*  expint  ( (b*g^2) . / (k (N) . ^2*U^4) )  only  for  PM  spec 


C.  MATLAB  CODE  FOR  APPLYING  SNELLS  LAW  TO  ELFOUHAILY 
SPECTRUM  MEAN-SQUARED  SLOPE 


%  09  sep  09 
%  apply  snells  law 
%parameters 

nl  =  1.3;  %  index  of  refraction  for  air 
n2  =  1 ;  %  index  of  refraction  for  water 

scale  =  0.01; 


R  =  y;  %  slope  (R  degrees  from  horizontal) 

thl  =pi/4;  %  angle  of  incidence 

th2  =  asin ( (nl /n2 ) *  sin (thl+R) ) -R;  %  angle  of  exiting 
Vth2  =  var(th2) 

VR  =  var(R) 

Y  =  Vth2/  VR 

%D  =  (1 . /sqrt (1-0 . 7^2  .*(sin(thl  +  R) ) . ^2 ) ) . *0 . 7 . *cos (thl  +  R)  -1 
%E  =  D.^2 
plot (thl , Y, ' X ' ) 
hold  on 

title ('ratio  of  variances  vs  thl  for  wind  speed  10') 
gtext('N  =  2^18,  Km  =  160') 
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